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f~i ' Abstract. We report on a Monte Carlo method that generates one-dimensional 

^pH, trajectories for Bohm's formulation of quantum mechanics that doesn't involve 

■y ' differentiation or integration of any equations of motion. At each time, t = nSt 

j^ ■ (n = 1, 2, 3 . . .), A^ particle positions are randomly sampled from the quantum 

^ I probability density. Trajectories are built from the sorted N sampled positions at 



each time. These trajectories become the exact Bohm solutions in the limits N ^ 00 
and St -^ 0. Higher dimensional problems can be solved by this method for separable 
wave functions. Several examples are given, including the two-slit experiment. 
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1. Introduction 

Given a solution to Schrodinger's wave equation, "ipix, t) = R{x, t)e^^^^'*'^/^ for real R{x^ t) 
and S{x, t), Bohm's formulation of quantum mechanics [H [21 [3] defines the equation of 
motion for particle trajectories in one-dimension as, 

dV dQ 

ox ox 
where V is the classical potential, and the quantum potential Q = —{h^/2mR)d'^R/dx'^. 
An ensemble of particle trajectories will maintain a density that is equal to the quantum 
mechanical probability density p{x,t) = \ilj{x,t)\'^ if the particle velocities are further 
constrained to obey the so called guidance law, 



. _ 1 dS 
m dx 



x=x(t) 

Even though the quantum potential does have some explanatory usefulness, typically 
if ip{x,t) is known, one only needs to solve the guidance law to obtain the particle 
trajectories. The guidance law, in general, is impossible to solve in closed form and 
must be solved numerically [H H] . 

In Brandt et al. E] it was shown that the Bohm trajectories can be defined by a 
concept they have named quantile motion (a purely non-physical probability concept). 
A unique trajectory is obtained by requiring that the cumulative probability function 
(CPF) for the quantum probability density p{x,t) = \ip{x,t)\'^ be constant along the 
trajectory, 

rxpit) 

CPF{x) = P = / p(a;,t) (is = constant, (3) 



where the trajectory is given by xp{t) = CPF^^{P = constant). Alternatively, 
one could think of the constant CPF restriction as total left (or right) probability 
conservation. Like Bohm's equations of motion, the quantile motion solutions can not 
be obtained in closed form (except for particular probability densities), and therefore, 
again, must be solved numerically. 

In this paper, we present a Monte Carlo method for a given one-dimensional 
probability density, that approximates Bohm trajectories without solving either of 
Bohm's equations of motion [Eqs. ([T], [2])], or the inverse of the CPF in (J3]). In other 
words, our method does not need to take a derivative or an integral at any point in 
the process. This fact demonstrates clear computational advantage over the the typical 
numerical techniques used to solve Eqs. ([2]) and ([3]). The method takes two parameters: 
the number of particles in the ensemble N, and the size of each discrete time step 6t. It 
is shown that the approximate trajectories become the Bohm trajectories exactly in the 
limits N —>■ 00 and 6t ^ 0. A possible extension of the method into higher dimensions 
is discussed, and it is found that the method can be used to find Bohm trajectories for 
those cases of a separable wave function. Finally, we demonstrate the method for the 
harmonic oscillator, the free particle, and the two-slit experiment. Also provided is a 
two-dimensional infinite square well example for a separable wave function. 
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Figure 1. The Density Sampling Method. At each time t — nSt [n — 1, 2, 3, . . .), N 
points are sampled from the probability density p(x, t) and sorted. Trajectories are 
constructed by joining the i-th sorted point from each time step. 



P(x) 



rejected 




Figure 2. The von Neumann Acceptance-Rejection Method. For a given density p{x), 
a set of x-points is generated by uniformly placing random dots on the graph. If a dot 
is under the density curve, that dot's x value is place in the set. 



2. The Density Sampling Method 

2.1. Description 

The method assumes that a probabihty density, p{x,t), is known and given. The 
evolution of the density is depicted by an ensemble of A^ particles, see figure [TJ The 
trajectories of the N particles are constructed from locations at times t = n6t {n = 
1, 2, 3, . . .) with step size 6t. At each time, the probability density is sampled to generate 
a set of A^ possible x-points. The x-points are sorted numerically. The i-th trajectory 
in the ensemble is built from the i-th x-point of the sorted A^ points at each time 
step. Though there are many ways to generate a set of A^ points sampled from a 
given distribution [3 [8], in the examples below (section Hj), we use the von Neumann 
acceptance-rejection method [9], see figure [21 with a uniform proposal distribution pu. 
The quality of the trajectories generated relies on how one chooses the number of 
particles in the ensemble A^, and on the size of each time step St. We place the following 
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restrictions for these two parameters, 

A^ > 2Lpi„ax and N6t ^ -, (4) 

where L is the width of the domain of the x coordinate (generally this is limited to a 
range of values where the density is essentially non-zero), e is a small number < e <^ 1 
with the dimensions of a speed, and Pmax is the maximum value of the density for all 
positions x E L and for all times between the initial and final times. 

2.2. Connection to Bohmian Mechanics 

The method constructs the i-th trajectory from the i-th N sampled and sorted points 
at each time step. The particle trajectories, therefore, do not intersect by design (a 
familiar property of Bohmian trajectories). Hence, between successive time steps the 
approximate size of the maximum change in position 6x is of the order 6x ~ 2L/N. 
Identifying the density sampled trajectory as XDs(t), and the Bohm trajectory as XBit), 
we assume at t = that xds{0) = xb{0). We now describe the density sampled 
trajectory as Xosif) = XB(t)+6x(t) for some function 6x(t) < 2L/N such that 6x{0) = 0. 
Computing the cumulative probability function (CPF) value ([3]) for the density sampled 
trajectory to first order in 6x, 

/XDs(t) 
p{x,t) dx ^ Pb + pixBit),t)5x{t). (5) 

-oo 

Recall, the CPF value for the Bohm trajectory Pb is constant. From the restrictions and 
assumptions above, Pmax^x -C 1, therefore, p6x < Pmax^x ^ 1, so the density sampled 
trajectory fluctuates about the Bohm trajectory. 

While the parameter N places the location of the particle close to the actual 
Bohmian location at each time, the other parameter 6t (the size of each time step) 
fixes the density sampled speed approximately equal to the Bohm speed. The difference 
in the two speeds is on the order of 6x/{26t) , where again 6x is the size of the fluctuation 
about the Bohmian trajectory. From (jlj), we find that 6x/{26t) ^ e ^ 1, or that the 
difference in the speeds in quite small. We note that the density sampled trajectory will 
become the Bohm trajectory for N -^ oo and 6t -^ 0. 

3. Extension to Higher Dimensions 

In general, the one-dimensional density sampling method described above can not be 
extended into higher dimensions. In higher dimensions there is no natural ordering 
to sort the N sampled points, and therefore, no way to identify the i-th position 
in the ensemble at each time step like in the one- dimensional case. However, the 
quantile motion concept ([3]) can be used independently on each coordinate in higher 



dimensions when the wave function is separable as shown in [Appendix A[ Since the one- 
dimensional density sampled trajectory fluctuates about the Bohm trajectory defined 
by (E]), the one- dimensional density sampling method can be used independently on 
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each coordinate to generate higher dimensional Bohm trajectories for those cases of a 
separable wave function. We also note that the method will generate higher dimensional 
Bohm trajectories for wave functions that are nearly separable [in], by applying the 
method independently on each coordinate of the separable part of the wave function. 

4. Examples 

In the examples below, the probability density is determined in the usual way p{x,t) = 
\ip{x,t)\'^. The wave function was chosen to be non-stationary so that dp/dt ^ 0. For 
comparison, the Bohmian trajectories are computed from ([2]) or its equivalent, 

h fip*dip/dx — ipdip* /dx\ 
2^i \ ^ ) ■ ^ ^ 

In all cases, the range L of possible position values x was limited to an area where the 
probability density was essentially non-zero. From (j3]), the number of particles in the 
ensemble (or the number of sampled points) A^ was approximately equal to 2Lpinax x 10^, 
while e was taken to have a numerical value of the order 10~^. 

4.1. Harmonic Oscillator 

The harmonic oscillator wave function was a superposition of the ground state and the 
first three odd excited states, 

ip{x,t) = — ^^e ^^ > , (7) 

2v/^ „4^3,, Vd^ 

where a = y^h/miu, En = hhj{n + 1/2), and iJ„ are the Hermite polynomials. 
Naturalized units were used so that ^ = 1, u; = 3, and m = 1. The range of time 
was t G [0,3], and the size of each time step was 5t = 0.1. The range of the possible 
positions was x G [—5,5] (an area where the density was essentially non-zero). The 
number of particles in the ensemble was N = 10^. Five of the resulting density sampled 
trajectories (-I-) are shown in figure [3] against the actual Bohm trajectories. Notice that 
the density sampled trajectories are able to depict the oscillatory behavior of the Bohm 
trajectories rather well. 

4.2. Free Particle 

The wave function for the free particle (assumed to be Gaussian initially with width a) 
was taken to be, 

(9 \ 1/4 -ax^ /[l+{2ihat/m)] 
- \ (8) 

Naturalized units were used so that h = 1, m = 1, and a = 7r/2. Again, t G [0, 3] with 
time steps of St = 0.15. The number of particles in the ensemble was A^ = 10^. Six of 
the resulting density sampled trajectories are plotted in figure H] superposed on top of 
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Figure 3. (color online). Harmonic oscillator with wave function as a superposition 
of ground and the first three odd excited states. The Bohmian trajectories are solid 
lines while the density sampled trajectories are plotted as plus (+) signs. The plot is 
in naturalized units. 




Figure 4. (color online). Free particle with an initial wave function of a Gaussian 
centered around zero. The Bohmian trajectories are solid lines, and the density 
sampled trajectories are plotted as plus (+) signs. The plot is in naturalized units. 
The ensemble of trajectories demonstrates the familiar spreading of the wave packet. 



their corresponding Bohm trajectory. The trajectories depict the famihar spreading of 
the wave packet. 



4-3. Two- Slit Experiment 

This two-sht example is from §5.1.2 in Holland [3j. At first this problem seems to 
be two-dimensional. However, the motion along the coordinate from the slits to the 
screen [x in figure [5] is assumed uniform, thus the probability density is effectively one- 
dimensional. To allow for easier computation the experimental values given in Holland's 
book were rescaled so that h = 1, m = 1, and the total time between the slits and screen 
was tmax = 100. The range of possible positions was y G [—129.668, +129.668]. The 
number of particles in the ensemble was A^ = 10^, and 6t = tmax/30. The resulting 
trajectories were then rescaled back to Holland's numbers for plotting. Thirty of the 
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Figure 5. (color online). Two-Slit Experiment as described in §5.1.2 of P.R. Holland, 
The Quantum Theory of Motion: An Account of the de Broglie-Bohm Causal 
Interpretation of Quantum Mechanics, (Cambridge University Press, New York, 
1993). The density sampled trajectories (+) are plotted superposed on the Bohmian 
trajectories (solid). The initial positions are assumed Gaussian in the slits (left side of 
figure), and the ensemble of trajectories makes the familiar bands of bright and dark 
on the screen (right side of the figure) . 



density sampled trajectories (+) are plotted against their Bohm counterpart in figure [5l 
The trajectories manifest the familiar bright and dark bands of the two-slit intensity 
pattern on the screen (left side of figure). Also, notice the size of the fluctuations about 
the Bohm trajectory in the different regions. In high density regions the method does 
better since 6x oc 1/A^ is smaller. But in low density regions (between the bright bands) 
Sx is larger due to less particles being there. 



4-4- Separable Wave Function in a Two-Dimensional Infinite Square Well 

In figure [6] is a comparison of the density sampled trajectories and their Bohm 
counterpart for the two-dimensional infinite square well. The separable wave function 
was assumed to be ip{x,y,t) = ipx{x,t)iljy{y,t), where, 



'4'xix,t) 



_.sin(^)e-^^i*/^ + sm 



27rx 



-iAEit/h 



(9) 



and a similar expression for ipy{y,t). The energy Ei = n h /2mL , and naturalized 
units were used so that m = 1, h= 1, and the width of the well in each direction taken 
to be L = 1. The one- dimensional density sampling method was used independently 
for each coordinate. The time was t = n6t G [0, 1] for n = 1, 2, 3 . . ., and the size of 
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Figure 6. (color online). Comparison of the density sampled (+) and the Bohm 
trajectories for a separable wave function in the two-dimensional infinite square well 
of size 1 X 1 in naturalized units. The initial position of each trajectory lies on the line 
from (0.5,0) to (0,0.5). 



each time step was 6t = 0.05. The number of particles in each coordinate's ensemble 
was A^ = 10^. Note, in general, a particle's place in the each coordinate's ensemble is 
not the same. The density sampled trajectory points (+) are plotted superposed on top 
of the corresponding Bohm trajectories. For the separable wave function, the density 
sampled trajectories are again identical to the Bohm trajectories. 



5. Concluding Remarks 

The density sampling method described in this paper shows that one can generate 
approximate Bohmian trajectories without differentiation or integration of any equations 
of motion. The method merely evaluates the given probability density at various 
locations and times. It uses two adjustable parameters. The first being the size of 
the set of sampled points A^, and the second the size of each time step 6t. It was shown 
that as A^ ^ cxD and 6t —* the density sampling method trajectories become the 
Bohm trajectories. Like the quantile motion idea, this method can generate the unique 
one-dimensional Bohmian trajectories without any additional physical assumptions to 
standard quantum mechanics. In addition, this method can be used to construct one- 
dimensional Bohmian trajectories for densities where the wave function is unknown or 
doesn't exist (i.e., non-quantum problems). Further, if a known higher dimensional 
wave function is separable (or nearly separable), then the method can be applied 
independently to each coordinate, thus constructing higher dimensional Bohmian 
trajectories as well. 
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Appendix A. Extension of Quantile Motion into Higher Dimensions 

The extension of the quantile motion into higher dimensions is also discussed in Brandt 
et al. [5]. They show that instead of the total left (or right) probability being conserved 
in one-dimension, that in higher dimensions the probability is conserved inside a volume 
that is enclosed by a surface of Bohmian trajectories. This property, however, is not 
unique to Bohmian mechanics. Any velocity field x will conserve the probability inside 
the volume in configuration space since [Ij (chapter four), 

dp „ . , dJ ^„ . / , , N 

-f- = -pV ■ X and ^ = +JV • x, (A.l) 

at at 

where p is the probability density and J is the Jacobian such that the volume changes 

as dV{t) = JdVo. The probability inside this evolving volume is, 

Pin= [ pdV{t)= f pJdVo. (A.2) 

Which implies that dP\n/dt = 0. Therefore, any velocity field x will conserve the total 
probability inside a volume that in enclosed by a surface of trajectories solved from x. 
This is in contrast to what is found in the one- dimensional case, where the velocity field 
that satisfies the total left (or right) probability conservation is unique. 

However, the quantile motion concept can be used to generate trajectories in higher 
dimensions if one uses the marginal distribution for each coordinate analogously to the 
cumulative probability distribution (CPF) in one-dimension. Suppose the system can 
be described by A^ Cartesian coordinates, then the corresponding definition of ([3]) is, 

/•Xi{t) 

Pi= Pi{x^,t)dxi for i= l,2,...,Ar, (A.3) 



where pi{xi,t) is the marginal distribution for the i-th coordinate. Again, we assume 
that the particles are conserved so, 

%.tli,,.)^<, (A.4, 

Partially integrating this continuity equation we have that, 

"Xiit) |-+oo g 

.. — dxi...dxN 

oo J —oo J —oo '^'^ 



N 



i=i 



Xi(t) 



d{pxj 



+ 2^ ■■ ■■ -^^dxi...dxN = 0. (A.5) 



"^3 



Interchanging the partial derivative with respect to time and performing the ±cxo 
integrations, and assuming that the density is zero at ±oo, only the j = i integrals 
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on the second term survive, 

"""' ^ dx 
, ot 

/ + 00 p+oo 

.. / pxi dxi . . .dx^i . . .dxN = 0. (A. 6) 

-oo J —oo 

Then with the expression above, and differentiating ( 1A.3P with respect to time. 



PiXi — / .. / pxi dxi . . .dx^i . . .dxN- (A. 7) 

'^'' J-oo ^— CO 

Hence, the i-th coordinate, in general, doesn't conserve total left probability of the 
marginal distribution since Xi could depend on the other coordinates. Suppose, however, 
that Xi = Xi{xi,t) (i.e. the motion along the i-th coordinate is independent), then 
dPi/dt = 0, and the total left probability is conserved for the marginal distribution p,. 
In higher dimensional Bohmian problems the guidance law ([2]) for the i-th 
coordinate becomes, 

1 dS{xi,X2,...,XN;t) 

Xi = . (A. 8) 

m oxi 

If the wave function is separable, then, 

ip = ipi{xi, t)tp2{x2, t) ■ ■ ■ iPn{xn, t). (A. 9) 

The probability density is also separable, p = pi{xi,t)p2{x2,t) . . . p^{xj\f,t), and the 
phase becomes S = Si{xi,t) + S2{x2,t) + . . . + S^i^x^jt). Which implies that 
Xi = {l/m)dSi{xi, t)/dxi, or that the velocity field for the i-th coordinate is independent 
of the other coordinates, for alH = 1, 2, . . . , A^. Hence, one can use the one-dimensional 
total left (or right) probability conservation, independently for each coordinate, to 
generate higher- dimensional Bohmian trajectories for a separable wave function. 
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